calcBPdirectory <- function(directory = "", comp.rands = 2, ...){
  
  ldir <- length(dir(directory))
  df.out <- data.frame(rand = NA, sim.sb = NA, sim.cr = NA, d2.p = NA, hz = NA, vdh = NA, rmse = NA)    
  
  for(i in 1:ldir){
    load(paste(directory, "sim", i, ".RData", sep=""))
    sb.data <- bdata$orig
    stats.sb <- calcBalancePrecision(sb.data, ...)
    
    df.out <- rbind(df.out, c("sb", i, NA, stats.sb$d2.p, stats.sb$hz, stats.sb$vdh, stats.sb$rmse))
    
    for(cr.idx in 1:comp.rands){
      
      sb.data$Tr <- sample(c("Treatment 1", "Treatment 2"), nrow(sb.data), replace = TRUE)
      stats.cr <- calcBalancePrecision(sb.data, ...)
      
      df.out <- rbind(df.out, c("cr", i, cr.idx, stats.cr$d2.p, stats.cr$hz, stats.cr$vdh, stats.cr$rmse))
    }
    
  }  
  
  ## remove leading NA row:
  df.out <- df.out[2:nrow(df.out), ]
  
  df.out <- reshape(df.out, varying = list(c("d2.p", "hz", "vdh", "rmse")) , v.names = "stat.val", direction = "long")
  df.out$stat.name <- recode(df.out$time, "1='d2.p'; 2='hz'; 3='vdh'; 4='rmse'")
  
  df.out$sim.sb <- as.numeric(df.out$sim.sb)
  df.out$sim.cr <- as.numeric(df.out$sim.cr)
  df.out$stat.val <- as.numeric(df.out$stat.val)
  row.names(df.out) <- 1:nrow(df.out)
  
  return(df.out)
}